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I . INTRODUCTION 


This final report on NASA Grant NGL 33-010-070, "Finite 
Element Shell Instability Analysis" summarizes the work accom- 
plished under the grant from its inception in November, 1968 until 
1 July, 1975, the date of termination of the grant. The overall 
objectives of this work were to formulate procedures and an asso- 
ciated computer program for finite element thin shell instability 
analysis* In order to meet these objectives, the research was 
divided into the following component aspects: 

(a) formulation of basic element relationships 

(b) construction of solution algorithms on both the 
conceptual and algorithmic levels 

(c) conduct of numerical analyses to verify the accuracy 
and efficiency of the theory and related programs. 

These aspects and the results achieved therein are described 
in subsequent sections of this report. First, however, a section 
is devoted to a brief description of the background of finite 
element thin shell instability analysis. The body of the report 
concludes with remarks concerning theoretical directions and 
numerical schemes, identified in the course of the grant research, 
which could not be pursued in detail in conjunction with the grant 
but which hold much promise for improvements in finite element 
thin shell instability analysis. 

NASA Contractor Reports produced under this grant, as well as 
more informal reports and papers published or presented with 
support of the grant (or deriving perspectives therefrom) are 
represented in the References of this report. 
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II. BACKGROUND 

The problem of finite element thin shell instability analysis 
represents the confluence of three major technological disciplines, 
each of which has been the subject of intensive and widespread 
study in past years. These disciplines are thin shell theory, 
structural instability theory, and finite element analysis. Other 
disciplines have been also involved, such as numerical methods 
(eigenvalue calculation, numerical integration) and computer 
programming, but this is usually regarded as either subsidiary or 
peripheral to those delineated above. 

In thin shell theory the clarification of appropriate strain- 
displacement relationships has been a relatively recent development 
(since 1960). Certain formulations have been known to be deficient 
in that they give rise to strain under rigid body motion. Other 
formulations do not possess such shortcomings and differ from one 
another only with respect to terms which are unimportant in com- 
parison to the assumptions made in formulating a thin shell theory. 
These have been identified in papers by Koiter (Ref. 1) and by 
Sanders (Ref. 2) and characterized as '‘consistent" formulations. 

The present work has basically adopted the consistent formulation 
due to Koiter (Ref. 1). Koiter's equations apply only to linear 
analysis; shell instability analysis requires the inclusion of 
terms that represent geometrically nonlinear behavior. The terms 
formulated by Mushtari and Galimov (Ref. 3) are appropriate to 
the latter and have been adopted in the work described herein. 

Structural instability theory has sustained a similar growth 
of understanding. The parameters to be solved for differ from one 
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physical circumstance to another and they do not pertain exclu- 
sively to a tracing of the load-displacement response along a 
continuous curve, starting from the linear regime at small loads 
into significant nonlinearities at higher loads. Rather, the 
diverse situations portrayed in Figure 1 are encountered. This 
figure examines the load-displacement behavior of different types 
of structural action, tracing in each case the response of a 
representative degree-of- freedom. 

The solid line of Figure 1 applies to "perfect" structures 
and represents the case in which the structure first displaces 
along the path defined by OAB (the fundamental path) and bifurcates 
(or branches ) at the point A. The postbuckling path may rise 
(path AC), as in the case of compressed perfectly- flat plates, or 
it may descend (path AD), as in the case of compressed, perfect, 
thin-walled cylinders with appropriate support conditions. 

When the structure possesses fabricational imperfections, the 
load-displacement behavior in the presence of destabilizing phe- 
nomena follows the paths indicated by dotted lines (Figures la and 
lb). Structures with a rising post-buckling path as in Figure la 
(e.g., flat plates) will have strength exceeding the bifurcation 
load. The strength of an imperfect structure with a descending 
post-buckling path in the perfect state, represented by 
Figure lb, will not achieve strengths as high as the bifurcation 
load unless the load-displacement path again rises at larger dis- 
placements. Such structures, under the appropriate load condition, 
are therefore imperfection sensitive and the maximum load attained 


(Point E) is termed the limit point . 
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A i.on*bifurcating load>displacement behavior may also occur 
for a structure assumed to be devoid of imperfections and may take 
the form similar in shape to the curve OE (Figure lb) of the 
imperfection-sensitive structure. Again, a limit point is encoun- 
tered and the buckling phenomenon is of the snap -through type. 

Four component areas of structural instability theory may 
therefore be perceived from the foregoing: (1) general nonlinear 

prebuckling analysis, (2) the calculation of limit points, (3) cal- 
culation of bifurcation (or ^branching") points, and (4) determina- 
tion of the load-displacement response along a post-buckling path. 

Surveys of the finite element method in structural instability 
theory have appeared at various times since 1969, when Martin 
(Ref. 4) presented a view of the topic. More recent surveys have 
appeared in Refs. 5-8. 

In regard to finite element analysis it is fair to say that 
no acceptable curved thin shell element was available at the start 
of the current research. The preferable formulation would permit 
accurate analyses at reasonable cost and be reliable in the analysis 
of the widest range of shell forms and loadings. The desirability 
of meeting all relevant conditions on a finite element formulation- - 
but especially those on zero- strain under rigid body motion, inter- 
element continuity of displacement, and the representation of 
constant strain- -is widely accepted. All of these cannot be 
achieved at reasonable computational cost, however, and the inter- 
vening years have seen extensive studies of the tradeoffs that can 
be made. A recent summary of the state-of-the-art in finite 
element thin shell analysis has been given in Ref. 9. 
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III. ELEMENT FORMULATION 

Major components of the efforts on this grant were devoted to 
the formulation of the finite elements needed for instability 
analysis of general curved shells. The basic element is an arbi- 
trary triangular shell element (Fig. 2), although formulations 
were accomplished for an arbitrary quadrilateral element (Fig. 3) 
and a stiffening element (Fig. 4). 

As noted in the prior section, the entire field of finite 
elements for shell analysis was in its infancy at the start of the 
work on this grant. Although linear shell analysis is well in 
hand, the amount of work done with respect to elastic instability 
analysis in this mode continues to be quite limited. 

The usual conditions on a finite element formulation (compati- 
bility, constant strain, rigid body motion) have been approached, 
if not met exactly, through formulation of the triangular element 
on the basis of the Koiter (Ref. 1) ’’consistent" shell theory and 
cubic polynomial descriptions of the three displacement components. 
This formulation does not have the required interelement continuity 
of displacement, so continuity is "restored" by writing constraint 
conditions which are appended to the analysis with use of the 
Lagrange multiplier method. 

The formulation of the above element, which is done in terms 
of orthogonal curvilinear coordinates and with nonzero Gaussian 
curvature, is described in detail for linear analysis in CR-2482 
(Ref. 10). This report includes extensive numerical solution data 
which demonstrates that the approximations made with respect to 
the condition on zero-strain under rigid body motion are satisfac- 
tory for all of the conventional problems studied as well as for 
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the particular pathological cases studied. The latter were intended 
to emphasize any deficiencies which might arise on account of the 
failure to meet exactly the condition on zero-strain under rigid 
body motion. 

The extension of the triangular thin shell element to account 
for geometrically nonlinear behavior is detailed in CR-2483 (Ref. 
11). Here, the shell equations are extended to the description of 
nonlinear behavior with terms derived by Mushtari and Galimov 
(Ref. 3). The finite element representation of these terms was 
accomplished with the same displacement assumptions as for the 
linear terms ("consistent” formulation), as well as with simpler 
assumptions ("inconsistent" formulation). Also, a combination of 
four triangles was made in formation of a quadrilateral element. 

It should be noted that evaluatica of the strain energy 
integrals that give the element stiffness coefficients was performed 
by means of numerical integration, using Gauss -Radau quadrature. 
Numerical studies of certain problems, particularly those featuring 
spherical caps, disclosed that the rapid variation of geometric 
parameters in the vicinity of the crown could not be adequately 
handled by linear variation of these parameters within the element. 
By implementation of a capability to evaluate the geometric parame- 
ters at each numerical integration point it was possible to obtain 
highly accurate results for the problems in question. This more 
recent work is described in Ref. 12. 
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IV. SOLUTION ALGORITHMS - CONCEPTUAL 

It vas noted in Section II of this report that there are 
four divisions of the general topic of shell instability analysis: 
(1) nonlinear prebuckling, (2) bifurcation, (3) limit point, and 
(4) postbuckling. 

In the present work, in nonlinear prebuckling analysis, the 
approach taken was the incremental -iterative method. The load 
path is first divided into a number of increments. The tangent 
stiffness is formed for a typical interval (see Fig. 5) and the 
linearized analysis is performed. Then, the calculated displace- 
ments are substituted into the equilibrium equations to yield a 
•’residual", which represents the error due to the linearization 
of the analysis. The residual is employed in an "initial-force" 
corrective analysis for the interval; this step does not require 
re-formation of the stiffness in the interval. After the initial 
force correction is made the analysis proceeds to the next interval. 

The calculation of bifurcation buckling, for most of the 
grant, was performed by interpolation of the determinant of the 
nonlinear prebuckling analysis stiffness equations through the 
zero value. This is based on the fundamental condition that at 
buckling the stiffness matrix is singular, i.e., has a zero value 
for its determinant. Subsequent experience showed this procedure 
to be unreliable, since the plot of determinant versus load might 
be erratic. Thus, in the latter part of this research the Sturm- 
soquence algorithm of Gupta (Ref. 15) was implemented and adopted. 
This algorithm is foolproof in the definition of the lowest eigen- 
value. 
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To trace the load-deflection curve up to and past the lim4>W • 
point special precautionary measures have to be taken. At the 
limit point the determinant of the stiffness matrix is zero and, 
consequently, the system of equations does not possess a unique 
solution. The deficiency is of rank one so to overcome this one 
constraint must be imposed on the equations. The present con- 
straint consists of the assignment of a prescribed value to one of 
the displacement degrees-of -freedom. 

The procedure is as follows. First, a representative degree- 
of-freedom is identified. The central displacement of a plate or 
shell is an appropriate choice, the selection being input directly 
by the user. Three load vectors are dealt with: 

(1) No loading, except for the reaction at the prescribed 

displacement: {Pj^} 

(2) The incremental load vector, with zero displacement 

at the representative freedom: {P 2 ) 

(3) The out of balance load vector calculated on the basis 

of the displacements at the last load level, with zero 
displacement at the representative freedom: {P^) 

These three load cases are shorn in Fig. 6 for an arch. The 
reactions at the constrained freedom are denoted by and R^, 

respectively. The desired combination of displacements is made 
up of two separate combinations thus 


and 



The first equation gives the necessary combination in the 
absence of out of balance forces. This would be sufficient in 
the case of the regular incremental approach. In that case 
would be zero at the limit point but R 2 would be nonzero and there 
would be no computational difficulty. 

In postbuckling analysis it is necessary to distinguish 
between limit point and bifurcation buckling. Limit point post- 
buckling can be handled by the procedures described above, i.e., 
by incrementing displacement rather than load. In bifurcational 
buckling it is necessary to extract the prebuckling path. This 
can be done by perturbation procedures, developed initially by 
Koiter (Ref. 14) for classical stability analysis and subsequently 
formulated by Thompson, et al (Ref. 15) for multiple degree-of- 
freedom (discrete coordinate systems) . A rather complete approach 
to finite element postbuckling analysis, based on perturbation 
procedures, was developed under this grant and is described in 
Refs. 16-18. 
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V. SOLUTION ALGORITHMS - NUMERICAL 

Efforts were devoted, during the full span of the grant, to 
the development of a computer program for finite element £hell 
instability analysis ("FESIA”). This program, coded in Fortran IV, 
is documented in Ref. 19, has been operational at both Cornell 
University and NASA Langley Research Center. Ref. 19 comprises 
a Level 3 program documentation as defined by NASA (Ref. 20). This 
means that theoretical background is outlined, operational, input, 
and output instructions are detailed and exemplified, and a flow 
chart of macro -operations is given. In view of the scope and 
content of Ref. 19 the nature of the program is only sketched in 
the following. 

The cornerstone of the program is the triangular thin shell 
element described in Section III. The element may have orthotrcpic 
material properties. The geometric description of the element 
permits non- zero Gaussian curvature referenced to a system of 
orthogonal curvilinear coordinates. Program options permit simpli- 
fied geometric description when common shapes (e.g., circular 
cylinders) are to be analyzed. The definition of support conditions 
is also simplified via the inclusion of familiar circumstances, 
such as fixed support. Point loads, distributed loads, and pre- 
assigned joint displacements can be specified. In the case of 
distributed loads a highly accurate description is possible, since 
the evaluation is performed at all numerical integration points 
within the element. It should be noted that a quadrilateral element 
can be specified which is formed of four triangular elements. 

The program permits a general, geometrically nonlinear analysis 
based on a Lagrangian description of strain-displacement equations 


i 


11 


and the incremental* iterative computational scheme. The load* 
displacement behavior can be carried beyond the limit point through 
use of displacement incrementation. Bifurcational buckling is 
determined through either an interpolation of the determinant or 
by use of Gupta's Sturm sequence algorithm (Ref. 13). 

A variety of analyses performed with use of the FESIA program, 
intended for comparison with alternative solutions, is described 
in the next two sections. 
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VI. VERIFYING ANALYSES 

In any large-scale computational approach it is essential that 
the adequacy of the formulation and of the underlying approximations 
be verified through performance of comparisons with alternative 
theoretical solutions and test data. Availability of such informa- 
tion is limited for the case of elastic instability of general 
curved thin shells. 

A rather extensive series of verifying analyses for the 
formulations developed under this grant has been reported in the 
relevant Contractor Reports and published papers. 

These include the "standard" comparison problems of the cylindrical 
shell roof originally analyzed by Scordelis and Lo (Ref. 21), the 
pinched cylinder (Ref. 22), a hyperbolic paraboloid shell (Ref. 23), 
a torus (Ref. 24), and a pinched and pressurized sphere (Refs. 25 
and 26). In all cases the solution accuracy and computational 
efficiency was found to be comparable, if not superior, to alter- 
native formulations. 

During the final year of the grant some comparison analyses 
were performed of two particular problems which have represented 
a long-standing challenge in the prediction of shell elastic 
instability. These were the torispherical shell under internal 
pressure and the hyperboloid under wind loading. 

In 1959 Galletly (Ref. 27) called attention to the fact that 
when torispheres are used as the closures of internally-pressurized 
tanks, circumferential compressive stresses are introduced which 
may lead to an elastic instability mode of failure. The problem 
is quite complex due to the rapidly-varying state of stress and 
the localized nature of the buckle modes. 
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Figure 7 shows the finite element idealization used in 
analysis of this problem. The capability of accurately represent- 
ing the rapidly changing geometric parameters near the crown was 
especially important. Figure 8 shows the level of agreement 
between an alternative solution for the prebuckling (linear) 
stresses (Ref. 24) and the finite element solution; these are seen 
to be in reasonable agreement. V^ith respect to buckling, the 
bifurcation solution obtained in this study was 6.68 p.s.i., 
compared with an experimental value of 6.7 in the experiment of 
Adachi and Benicek (Ref. 28). 

The torispherical head represented an especially difficult 
bifurcational analysis condition. The lowest positive eigenvalue 
corresponded to buckling under external pressure, a negative 
value. Also, the relationship between the load intensity and 
determinant was quite erratic, requiring use of the Sturm sequence 
algorithm (Ref. 13). 

The hyperboloid under wind loading (Fig. 9) also represents 
a problem with rapidly varying (spatially) prebuckling stresses. 
Much attention has been given to this problem because of the 
failure of a number of cooling towers under wind load in England 
in the mid-1960*s. It should be noted that a similar structural 
form is encountered in many aerospace applications (e.g., rocket 
engine thrust chambers) . 

A finite element idealization of the cooling tower is seen 
in Figure 9. Figure 10 gives a comparison of results obtained 
with the FESIA program (Ref. 12) and those due to Fonder (Ref. 26) 
and Cole, et al (Ref. 34). They are seen to be in close agreement 
with the 1 letter. A comparison of buckling test data (Ref. 30) 
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and the present solution showed a difference of 251 (Ref* 12). 
This is closer SL^reement than recorded in studies performed else- 
where. 

VII. THEORY- TEST COMPARISON 

One theory- test comparison accomplished under this grant, 
results of which have design significance, involved the shear 
buckling of reinforced, perforated plates (Figs. 11, 12). A 
practical question which has • ^t-^d for many years concerns the 

amount of reinforcing material that is needed to compensate for 
the material removed in making the hole. A number of plates of 
different proportions of reinforcement were tested under support 
of the Cornell Structural Engineering Department and analyses of 
these conditions were performed by means of the FESIA program. 

The theory-test comparisons were good and disclosed that the 
amount of reinforcement needed to compensate for the hole is con- 
siderably less than had heretofore been assumed. This work is 
described in a NASA CR (Ref. 31). 
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VIII. CONCLUDING REMARKS 

Not all of the theoretical work performed was implemented in 
the FESIA program as of the conclusion of the grant, and various 
schemes have been identified which hold much promise for the 
improvement of the accuracy, scope and efficiency of the program. 
These include: 

(Ij Treatment of interelement displacement discontinuity, 
in the case of the triangular element, can be accom- 
plished in a manner which eliminates the need for 
Lagrange multipliers. One such approach has been 
proposed by Kikuchi and Ando (Ref. 32). 

(2) Self-contained assessment of whether a problem is 
in the class of bifurcation buckling or limit point 
and succeeding computation strategy based on the 
result. Necessary theoretical basis is outlined in 
Ref. 33, a paper written under support of this grant. 

(3) Implementation of various ’’condensation” schemes, 
enabling global analysis for bifurcation and limit 
point behavior to be conducted for a reduced number 
of degrees-of-freedom. 

(4) Expanded representation of various types of elements, 
with special attention given to stiffeners and the 
verification of stiffened shell analysis, 

(5) Improved input-output facilities. Emphasis should 
be placed on computer graphics, with use of picture 
tube representations and including mesh generation 
capabilities. 
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Work on all of the above considerations, unsupported by 
grant funds, was in progress at the time of writing of this report 
and will be made available to NASA when complete, if so desired. 
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Figure 2. Triangular Thin Shell Element 
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Figure 3. Quadrilateral Element 
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Figure 7. Finite element idealization of torispherical 
head pressure vessel and junction between 
component shells 
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Figure 8. Normal displacements w for pressurized torispherical 
head vessel 
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Figure 9. 
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Present study 4x5 mesh 
V Fonder (Ref. 26) 20 x 1 mesh 
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Figure 10. Hdrlzontal de forma tj. on of cooling: tower 
under wind load at <^=0.0 and 0.240 



Figure 12. Reinforcement Patterns Around Perforations in 
Plates with d/b = 0.4 














